library(tidyverse)
library(PRISMAstatement)
Global urbanization continues unabated, with more than 50% of the worlds’ population living in cities. Cities are conventionally viewed as a threat to local biodiversity because natural habitat is replaced with development. However, more recently, there is greater acknowledgement from the public and private sectors that supporting local environments sustains critical ecosystem services, which in turn improves human health and biodiversity conservation. Consequently, urban planning and design has shifted towards green infrastructure (GI), such as green roofs and retention ponds, to increase connections between city and nature in an era of climate change. The contribution of GI to some ecosystem services has been proven (e.g. stormwater management, building cooling), but the contribution to biodiversity conservation remains unspecified. Using a systematic literature review, this project will (i) determine effect estimates that relate different GI types and characteristics to the impacts on natural systems, and (ii) compile relevant data to develop different implementation scenarios GI for Toronto and region. This study will inform natural system planning and improve quantification of GI on urban ecosystems. Findings from this research will have global ramifications that allow city planners to optimize GI implementation for sustainable development and decrease the impacts of cities on natural systems.
| date | task |
|---|---|
| June 18 | Begin meeting with staff and MacIvor lab to set out workplan |
| June 25 | Begin literature review and data extraction |
| July 2 | Aggregate available data for GI analysis in Toronto |
| July 3 | Complete meetings with TRCA staff on relevant considerations for the project |
| July 9 | Determine important parameters for modelling GI in Toronto |
| August 20 | Complete collection and review of relevant articles |
| September 3 | Conduct meta-analysis on available data |
| September 10 | Propose candidate models for quantifying GI effects for natural systems |
| Sept 24 | Model validation and begin writting manuscript |
| Oct 15 | Complete a draft of manuscript and finalize model |
A systematic literature search will be conducted using Web of Scienc for all articles between 1980 and 2018. This time frame was chosen because it captures the majority of the literature on green infrastructure and with articles that have some measured estimates. The review will include all studies globally. The intended purpose of this search is to capture all articles that have documented both green infrastructure implementation and a measure of natural systems. The search terms that will be used are:
*green infrastructure OR low*impact development OR Sustainable Drainage System* OR Water Sensitive Urban Design OR green*roof AND *diversity OR species OR ecosystem OR ecology OR habitat* OR co-benefit
These terms have returned 871 results (as of June 2018)
search1.1 <- read.csv("data/WOS-lit.csv")
search1.2 <- read.csv("data/WoSPart3-July_4_2018.csv")
net.difference <- anti_join(search1.2, search1.1, by = "DOI")
net.difference <- net.difference %>% select(Title, DOI) #to simplify for a look
nrow(net.difference) #count of number of differences from consecutive search
## [1] 182
## 182 papers to be added by including revised terms
## Select those articles and join with other dataset
net.difference <- anti_join(search1.2, search1.1, by = "DOI")
updated.search <- rbind(search1.1, net.difference)
#write.csv(updated.search, "data/WOS-lit.updated.csv")
Adding revised terms from July 3rd meeting added 182 papers Total articles returned = 1,053 (as of July 2018)
## Adding terms for naturalized pond and pollinator garden
search1.2 <- read.csv("data/WOS-lit.updated.csv")
search1.3 <- read.csv("data/WoSPart4-July_11_2018.csv")
net.difference <- anti_join(search1.3, search1.2, by = "DOI")
net.difference <- net.difference %>% select(Title, DOI) #to simplify for a look
nrow(net.difference) #count of number of differences from consecutive search
## [1] 0
## 213 papers to be added by including revised terms
## Select those articles and join with other dataset
net.difference <- anti_join(search1.3, search1.2, by = "DOI")
updated.search <- rbind(search1.2, net.difference)
#write.csv(updated.search, "data/WOS-lit.updated.csv")
Adding revised terms from July 3rd meeting added 213 papers Total articles returned = 1,224 (as of July 2018)
This steps includes a. checking for duplicating, b. reviewing each instance for relevancy, c. consistently identifying and documenting exclusion criteria. Outcomes include a list of publications to be used for synthesis, a library of pdfs, and a PRISMA report to ensure the worflow is transparent and reproducible. Papers were excluded with the following characteristics:
evidence <- read.csv("data/evidence.csv")
excludes <- evidence %>% group_by(reason) %>% count(exclude) %>% filter(reason!="")
ggplot(excludes, aes(x=reason, y=n)) + geom_bar(stat="identity") + coord_flip()
## frequency of study
year.rate <- evidence %>% group_by(Year) %>% summarize(n=length(Year))
## Completed so far
prog <- sum(evidence$exclude!="")
prog
## [1] 1266
## Remaining
total <- nrow(evidence)
total
## [1] 1266
setTxtProgressBar(txtProgressBar(0,total, style = 3), prog)
##
|
| | 0%
|
|=================================================================| 100%
Initial pass for relevant papers complete.
GI.type <- evidence %>% group_by(GI.type) %>% count(exclude) %>% filter(GI.type!="")
ggplot(GI.type, aes(x=GI.type, y=n)) + geom_bar(stat="identity") + coord_flip()
Representations of relevant GI types found in papers
Prisma report
## total number of papers found
nrow(evidence)
## [1] 1266
## number of papers found outside of WoS
other <- read.csv("data/other.sources.csv")
nrow(other)
## [1] 28
## number of articles excluded
excludes <- evidence %>% filter(exclude=="y")
nrow(excludes)
## [1] 1123
## relevant papers
review <- evidence %>% filter(exclude!="y")
nrow(review)
## [1] 143
## papers for meta
meta <- evidence %>% filter(meta.=="yes")
nrow(meta)
## [1] 117
prisma(found = 1255,
found_other = 28,
no_dupes = 1283,
screened = 1283,
screen_exclusions = 1140,
full_text = 143,
full_text_exclusions = 0,
qualitative = 143,
quantitative = 66,
width = 800, height = 800)
## Loading required namespace: DiagrammeR
The research questions we are exploring:
require(ggmap)
### Start with base map of world
mp <- NULL
mapWorld <- borders("world", colour="gray50", fill="gray50") # create a layer of borders
mp <- ggplot() + mapWorld
## colorblind-friendly palette
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7","#000000")
meta <- read.csv("data//evidence.csv")
meta <- subset(meta, GI.type!="")
## plot points on top
mp <- mp+ geom_point(data=meta , aes(x=lon, y=lat, color=GI.type), size=2) + scale_colour_manual(values=cbPalette)+
theme(legend.position="bottom", text = element_text(size=20))
mp
## Number of studies extracted from online data
occurdat<-list.files("data//MS.data",pattern=".csv$",full=T)
length(occurdat)
## [1] 76
## 70 Studies found with usable data for synthesis
## load master datasets
meta <- read.csv("data//Master.GI.Datasets.csv")
## Omit repo 3 because duplicated with study 1305 and remove repo-9 because not equivalent GI comparisons. Removed Repo 3 because compare roof with ground
meta <- meta %>% filter(Study != "repo-3" & Study!="repo-9" & Study!="repo-1")
## Drop relative abundance because difference = 0
meta <- meta %>% filter(Estimate!="Relative.Abundance")
## Load packages and functions
library(reshape2)
library(metafor)
source("meta.eval.r") ## Multiple aggregate
## Create Unique identifier column
meta2 <- meta
meta2[,"UniqueSite"] <- paste(meta2$Study, meta2$Taxa.simplified, meta2$GI.compare, meta$Estimate, sep="-")
meta3 <- meta2 %>% filter(Infrastructure != "natural") %>% filter()
## Use function to extract summary statistics for comparisons
## meta.eval arguments are (meta.data, compare, ids , stats)
Infra.compare <- meta.eval(meta3, Infrastructure, UniqueSite, Stat )
## Combine the lists into same dataframe
## Rename Columns in second dataframe
Infra.stat <- Infra.compare[[2]] ## extracted statistics
names(Infra.stat) <- c("UniqueSite","green_mean","green_sd","grey_mean","grey_sd","green_n","grey_n") ## rename columns to match
Infra.raw <- Infra.compare[[1]] ## calculated statistics from raw values
## Join two dataframes
meta.stat <- rbind(Infra.raw, Infra.stat[, names(Infra.raw)])
meta.ready <- escalc(n1i = grey_n, n2i = green_n, m1i = grey_mean, m2i = green_mean, sd1i = grey_sd, sd2i = green_sd, data = meta.stat, measure = "SMD", append = TRUE)
taxa.test <- c(rep("invertebrate",5),"amphibian","birds.bats","birds.bats","invertebrate","invertebrate","birds.bats","birds.bats","invertebrate")
GI.test <- c("wall","wall","roof","roof","pond","pond","wall",rep("roof",5),"roadside")
#random-effects meta-analysis for green infrastructure vs grey
m1 <- rma(yi=yi, vi=vi, mods=~factor(taxa.test), data = meta.ready)
summary(m1)
##
## Mixed-Effects Model (k = 13; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -26.8388 53.6776 61.6776 62.8879 69.6776
##
## tau^2 (estimated amount of residual heterogeneity): 10.9891 (SE = 5.0836)
## tau (square root of estimated tau^2 value): 3.3150
## I^2 (residual heterogeneity / unaccounted variability): 98.11%
## H^2 (unaccounted variability / sampling variability): 52.99
## R^2 (amount of heterogeneity accounted for): 11.12%
##
## Test for Residual Heterogeneity:
## QE(df = 10) = 128.5068, p-val < .0001
##
## Test of Moderators (coefficient(s) 2:3):
## QM(df = 2) = 3.8058, p-val = 0.1491
##
## Model Results:
##
## estimate se zval pval ci.lb
## intrcpt -1.7060 3.3264 -0.5129 0.6080 -8.2256
## factor(taxa.test)birds.bats -3.3826 3.7411 -0.9042 0.3659 -10.7150
## factor(taxa.test)invertebrate 0.6600 3.5310 0.1869 0.8517 -6.2607
## ci.ub
## intrcpt 4.8136
## factor(taxa.test)birds.bats 3.9498
## factor(taxa.test)invertebrate 7.5806
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Produce a forest plot to determine the effect sizes for each study
forest(m1, slab = meta.stat$UniqueSite)
## Check for publication bias
## The symetrical distriubtion suggests there is no publication bias
funnel(m1)
## Create Unique identifier column
meta2 <- meta
meta2[,"UniqueSite"] <- paste(meta2$Study, meta2$Taxa.simplified, meta2$Nat.compare, meta2$Estimate, sep="-")
## Remove comparisons except urban and rural
meta2 <- meta2 %>% filter(Habitat == "urban" | Habitat == "natural")
## Determine the number of comparisons available
compare.eval(meta2, Habitat, UniqueSite)
## [1] 20
## Use function to extract summary statistics for comparisons
## meta.eval arguments are (meta.data, compare, ids , stats)
nat.compare <- meta.eval(meta2, Habitat, UniqueSite, Stat )
## Combine the lists into same dataframe
## Rename Columns in second dataframe
nat.stat <- nat.compare[[2]] ## extracted statistics
names(nat.stat) <- c("UniqueSite","natural_mean","natural_sd","urban_mean","urban_sd","natural_n","urban_n") ## rename columns to match
nat.raw <- nat.compare[[1]] ## calculated statistics from raw values
## Join two dataframes
meta.stat <- rbind(nat.raw, nat.stat[, names(nat.raw)])
meta.ready <- escalc(n1i = urban_n, n2i = natural_n, m1i = urban_mean, m2i = natural_mean, sd1i = urban_sd, sd2i = natural_sd, data = meta.stat, measure = "SMD", append = TRUE)
#random-effects meta-analysis for urban GI vs natural
m1 <- rma(yi, vi, data = meta.ready)
summary(m1)
##
## Random-Effects Model (k = 20; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -14.8475 29.6950 33.6950 35.5839 34.4450
##
## tau^2 (estimated amount of total heterogeneity): 0.0758 (SE = 0.0539)
## tau (square root of estimated tau^2 value): 0.2753
## I^2 (total heterogeneity / total variability): 49.98%
## H^2 (total variability / sampling variability): 2.00
##
## Test for Heterogeneity:
## Q(df = 19) = 39.3013, p-val = 0.0040
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.0889 0.0957 -0.9295 0.3526 -0.2764 0.0986
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Natural vs Urban GI
## Produce a forest plot to determine the effect sizes for each study
forest(m1, slab = meta.stat$UniqueSite)
## Check for publication bias
## The symetrical distriubtion suggests there is no publication bias
funnel(m1)
## The area of green infrastructure
names(meta)[20:21] <- c("GI.area","height")
meta.area <- subset(meta, GI.area>0)
## Determine unique identifier
meta.area[,"UniqueSite"] <- paste(meta.area$Study, meta.area$Taxa.simplified, meta.area$Nat.compare, meta.area$Estimate, sep="-")
## Summarize average richness by area sizes
area.stat <- meta.area %>% filter(Stat=="count" | Stat=="mean") %>% filter(Estimate=="richness") %>% group_by(GI.area, GI.type) %>% summarize(val=mean(Value)) %>% data.frame(.)
area.rich <- area.stat %>% filter(GI.type=="green roof" | GI.type=="green wall" | GI.type=="retention pond" | GI.type=="yards/home gardens" | GI.type=="public/community gardens")
library(ggplot2)
## Species richness per area
ggplot(area.rich, aes(x=GI.area, y=val, color=GI.type) ) + geom_point(size=3) + theme_bw() + scale_color_brewer(palette="Set2") + ylab("Average Species Richness") + xlab("Average Area of Green Infrastructure")+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
## Summarize average abundance by area sizes
area.stat <- meta.area %>% filter(Stat=="count" | Stat=="mean") %>% filter(Estimate=="abundance") %>% group_by(GI.area, GI.type) %>% summarize(val=mean(Value)) %>% data.frame(.)
area.abd <- area.stat %>% filter(GI.type=="green roof" | GI.type=="green wall" | GI.type=="retention pond" | GI.type=="yards/home gardens" | GI.type=="public/community gardens") %>%
filter(GI.area<50000) ## keep numbers approximately similar - removed outlier of 200,000
## Species richness per area
ggplot(area.abd, aes(x=GI.area, y=val, color=GI.type) ) + geom_point(size=3) + theme_bw() + scale_color_brewer(palette="Set2") + ylab("Average Abundance") + xlab("Average Area of Green Infrastructure")+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
## Compare pH for abundance
ph.data <- meta %>% filter(pH>0) %>% filter(GI.type=="retention pond" | GI.type=="natural water")
ph.abd <- ph.data %>% filter(Estimate=="abundance" & Stat=="count")%>% group_by(pH, GI.type) %>% summarize(val=mean(Value)) %>% data.frame(.)
## plot richness against pH
ggplot(ph.abd, aes(x=pH, y=val, color=GI.type) ) + geom_point(size=3) + theme_bw() + scale_color_brewer(palette="Set2") + ylab("Average Richness") + xlab("Average pH of Green Infrastructure")+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
## Compare pH for occurrence
ph.occ <- ph.data %>% filter(Estimate=="occurrence")%>% group_by(pH, GI.type) %>% summarize(val=mean(Value)) %>% data.frame(.)
m1 <- glm(Value~ pH, family=binomial, data=ph.data %>% filter(Estimate=="occurrence"))
summary(m1)
##
## Call:
## glm(formula = Value ~ pH, family = binomial, data = ph.data %>%
## filter(Estimate == "occurrence"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.8915 -0.8488 -0.8377 1.5229 1.6131
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.22413 1.45261 -0.154 0.877
## pH -0.08643 0.20324 -0.425 0.671
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 308.55 on 251 degrees of freedom
## Residual deviance: 308.37 on 250 degrees of freedom
## AIC: 312.37
##
## Number of Fisher Scoring iterations: 4
## plot richness against pH
ggplot(ph.occ, aes(x=pH, y=val, color=GI.type) ) + geom_point(size=3) + theme_bw() + scale_color_brewer(palette="Set2") + ylab("Average Occurrence") + xlab("Average pH of Green Infrastructure")+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
## plot richness against height
h.rich <- meta %>% filter(height>0) %>% filter(Estimate=="richness") %>% filter(Stat=="count" | Stat=="mean") %>% group_by(height, GI.type) %>% summarize(val=mean(Value)) %>% data.frame(.)
ggplot(h.rich, aes(x=height, y=val, color=GI.type) ) + geom_point(size=3) + theme_bw() + scale_color_brewer(palette="Set2") + ylab("Average Richness") + xlab("Average Height of Green Infrastructure")+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
## plot abundance against height
h.abd <- meta %>% filter(height>0) %>% filter(Estimate=="abundance") %>% filter(Stat=="count" | Stat=="mean") %>% group_by(height, GI.type) %>% summarize(val=mean(Value)) %>% data.frame(.)
ggplot(h.abd, aes(x=height, y=val, color=GI.type) ) + geom_point(size=3) + theme_bw() + scale_color_brewer(palette="Set2") + ylab("Average Abundance") + xlab("Average Height of Green Infrastructure")+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
## Pond salinity
salt <- subset(meta, Salinity>0)
length(unique(salt$Salinity)) ## not enough samples for a meaningful comparison
## [1] 4
## Pond depth
deep <- subset(meta, depth..m.>0)
length(unique(deep$depth..m.)) ## not enough samples for a meaningful comparison
## [1] 4
gf.data <- read.csv("data//GI.data//GreenRoofGeocoded.csv")
garden.data <- read.csv("data//GI.data//ComGardensdata.csv")
pond.data <- read.csv("data//GI.data//RetentionPondsGPS.csv")
## Extract coordinates only for GI
gf.data <- gf.data[,c("lon","lat")]
garden.data <- garden.data[,c("lon","lat")]
pond.data <- pond.data[,c("lon","lat")]
## combine into single dataset
GI.data <- rbind(gf.data,garden.data, pond.data)
## add column for data type
GI.data[,"GI.type"] <- c(rep("Green Roof",nrow(gf.data)),rep("Community Garden",nrow(garden.data)),rep("Retention Pond", nrow(pond.data)))
## interactive map
library(leaflet)
## Warning: package 'leaflet' was built under R version 3.4.4
m <- leaflet(data=GI.data) %>%
addTiles() %>%
addCircleMarkers(~lon, ~lat, color = ifelse(GI.data$GI.type == "Retention Pond", 'blue', 'red'), radius=5)
m